Anomalous finite-size effect in superconducting Josephson junction arrays 
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We report large-scale simulations of the resistively-shunted Josephson junction array in strip 
geometry. As the strip width increases, the voltage first decreases following the dynamic scaling 
ansatz proposed by Minnhagen et al. [Phys. Rev. Lett. 74, 3672 (1995)], and then rises towards 
the asymptotic value predicted by Ambegaokar et al. [Phys. Rev. Lett. 40, 783 (1978)]. The 
nonmonotonic size-dependence is attributed to shortened life time of free vortices in narrow strips, 
and points to the danger of single-scale analysis applied to a charge-neutral superfiuid state. 
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Dissipation of the supercurrent through superfiuid and 
superconducting films is caused by the flow of free vor- 
tices transverse to the current |l| || . This phenomenon 
was analyzed in detail by Ambegaokar, Halperin, Nel- 
son, and Siggia (AHNS) B. In the context of a uniform 
superconducting film in zero magnetic field, their theory 
predicts a power-law current-voltage (/-V) relationship 
with a temperature-dependent exponent. Despite the 
relatively simple construct of the theory, its experimen- 
tal and numerical verification has remained controversial 
|p[-pJ3|. It has been noted that boundary and finite-size 
effects can dominate the measured voltage drop across 
the sample at sufficiently low temperatures, making it 
difficult to achieve an unambiguous comparison between 
theory and experiments. 

An alternative theory of vortex flow dissipation which 
also yields nonlinear I-V curves in the superconducting 
state derives from the dynamic scaling hypothesis fn].[L2| . 
For superconducting films and Josephson junction arrays 
(JJA), the two approaches differ in their predictions of 
the exponent a(T) characterizing the power-law behavior 
V ~ I a below the Kosterlitz-Thouless-Berezinskii (KTB) 
transition temperature Tkt- According to the AHNS 
theory, 



aAHNsCT) = X + L 



(1) 



set of vortex-antivortex pairs, thereby invalidating the 
AHNS treatment. Nevertheless, the AHNS theory should 
still apply at sufficiently weak currents. This picture has 
not been borne out by a recent numerical study |l0| which 
shows persistently larger value of a than Eq. |l]) predicts. 

In this Letter, we present I-V results from large-scale 
simulations of the resistively shunted Josephson-junction 
(RSJ) array. A rectangular strip geometry is used to de- 
couple boundary effects introduced by the current leads 
from finite-size effects arising from vortex/antivortex mo- 
tion transverse to the current /. The main finding of 
our work is the existence of three distinct regimes as the 
strip width L y is varied. For L y < ~ I , one is 
in the linear-response regime where the near-equilibrium 
dynamic scaling hypothesis applies. In this regime, the 
voltage V decreases as L y increases, and reaches a min- 
imum value predicted by Minnhagen et al. When L y 
exceeds a second length scale Z r ~ I~ x l 2 , V saturates 
to the asymptotic value predicted by the AHNS theory. 
The intermediate regime lh < L y < 1? is characterized by 
a rising V as L y increases. This behavior is explained in 
terms of "self-recombination" of vortex-antivortex pairs 
under periodic boundary conditions. 

We begin with the model Hamiltonian of a two- 
dimensional (2D) JJA in zero magnetic field, 



where \ — t^JrI^bT and Jr is the renormalized spin- 
wave stiffness. On the other hand, the dynamic scaling 
analysis of Minnhagen et al. |0] suggests 
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(3) 



aMWJoCH = 2 X - 1. 



(2) 



The two expressions coincide at Tkt where x = 2, but 
the difference grows rapidly as one moves to low temper- 
atures. While there seem to be ample numerical support 
for the Minnhagen et al. scaling |l(|[l2]-[l4| , agreement 
with the AHNS theory has also been reported [§J. To 
rationalize the two scenarios, Bormann p5[ made an in- 
teresting suggestion that strong current creates a dense 



Here <pi denotes the phase of the superconducting or- 
der parameter on grain i, and J sets the strength of the 
Josephson coupling between neighbouring superconduct- 
ing grains. The summation in Eq. ([}]) is over all nearest 
neighbor grain pairs. In the RSJ dynamics, the network 
is insulated from the substrate, but normal currents are 
allowed between neighboring grains, in addition to the 
supercurrent. The dynamics of the phase variables fol- 
lows the Josephson relations plus the Kirchoff law for 
current conservation at each node |j, 
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Here u is the conductivity of normal links, and I Gy ± is 
the external current source which is nonvanishing only 
at the boundaries of the network. To model the effect 
of temperature T, a gaussian noise current r)ij is added 
to each link between two adjacent nodes i and j, with 
(%(*)»7i'i'(*')> = 2ak B T5i jt i> r d(t - t'). 

To integrate the coupled equations (|j), one needs to 
first solve for Vi = (h/2e)d4>i/dt. This can be done effi- 
ciently using a pseudo-spectral algorithm. In our simula- 
tions, the network is chosen to be a strip of L x columns, 
with L y nodes in each column. Open boundary condi- 
tions are used in the x-direction (along the strip) to allow 
for current inputs and outputs. Periodic boundary con- 
ditions (PBC) are used in the y-direction perpendicular 
to the strip. After performing the Fourier transform in 
the y-direction, the left-hand-side of Eq. (||) reduces to a 
tridiagonal form which can be easily solved by the Gauss 
elimination method. The time stepping is done using 
a second-order Runge-Kutta method with At = 0.1 (in 
units of ah/2eI Cl where I c — 2eJ/h is the critical current 
through each link) To minimize the boundary ef- 

fects, voltage drop across the sample is measured in the 
interior of the network. For the data presented below, 
we typically choose L x = 2048 + 2 x 128 and discard 128 
columns at each end for this purpose. 

When the external current I into each boundary node 
is set to zero, we recover a dynamical version of the XY 
model whose equilibrium properties have been well char- 
acterized. In particular, the KTB transition tempera- 
ture has been determined to be Tkt — 0.89 J/k B [ p"7[ . 
We have confirmed this value in the simulation by ex- 
amining the correlation function of the superconducting 
order parameter and the helicity modulus (or Jr), which 
exhibit universal behavior at the transition. Direct mea- 
surement of the I-V curve at this temperature has also 
confirmed the predicted exponent a (Tkt) = 3. 

We now describe the main results of our simulation. 
Figures 1(a) and 1(b) present the I-V data at T = 
0.8J/kB and T = 0.7 J/k B , respectively, for four different 
array sizes: L y — 8, 32, 128, and 512. Here the current I 
is measured in units of I c , and V, the voltage drop per 
column, in units of h/2e. As the current / decreases, the 
finite size effect becomes more and more evident. Exam- 
ining the data at the smallest / shown, we see a clear 
nonmonotonic size dependence: for small values of L y , V 
decreases with increasing L y . However, this trend is re- 
versed as L y is increased further. The solid and dashed 
lines in the figure correspond to the power-law behav- 
ior predicted by AHNS [Eq. (g)] and Minnhagen et al. 
[Eq. (g)], respectively, using X = 2.8 at T = 0.8J/k B 
and x = 3.5 at T = 0.7 J / k B . These values of x are ob- 
tained from measurements of the helicity modulus at the 
respective temperatures, and agree with previous stud- 



ies §. For T — 0.8 J/k B , the I-V data at L y = 512 is 
well-described by the AHNS theory. On the other hand, 
the envelop of data sets at smaller values of L y follows 
the Minnhagen et al. formula (g). For T = 0.7J/k B , 
we are only able to obtain reliable data for I > 0.1/ c 
due to the rapid decrease of V with decreasing /. Nev- 
ertheless, better agreement with the AHNS prediction 
is apparent on the low current side for the L y = 512 
array. The larger apparent slope of the data at higher 
currents may be attributed to the larger value of J e g due 
to the smaller length scale probed, which yields a higher 
Xcff = TfJ c s/k B T. 



10" 



10" 



10- 



> 10"' 



10- 



10"' 



10"' 



; (a)r=0.8/A B 


AD 




A/ / 




A8 






□ 32 




\ a/ 


128 
• 512 









(b) T=0JM B 




10" 



10" 



///. 





□ 


A 8 


• 




□ 32 


O' 




128 


□ 




• 512 








10"' 


/// 

c 





10" 



FIG. 1. Current- voltage curves for arrays of various strip 
width L y at two temperatures below Tkt: (a) T — 0.8J/fcs, 
and (b) T = 0.7 J/ ks- The value of L y for each data 
set is given in the legend box. Also shown are the pre- 
dicted power-laws by AHNS (solid line) and by Minnhagen 
et al. (dashed line) at each temperature. 
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FIG. 2. Nonmonotonic finite-size dependence of the volt- 
age on the strip width for three different sets of current and 
temperature values. 
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Figure 2 shows three sets of voltage data against trans- 
verse array size L y . Data at T = 0.8J/fcs are represented 
by open (/ = 0.1/ c ) and solid circles (I = 0.05/ c ), while 
the set at T = 0.7 J/ ks and I — 0.1I C are presented 
by open squares. Statistical error of each data point is 
smaller or comparable to the symbol size. As / decreases, 
the size l m where the minimum value V — V m is reached 
increases. l m also varies with T, though this dependence 
may be weak. Since the change of V with L y is slow in 
the neighborhood of l m , it may give the false impression 
that the finite-size effect is already saturated at l m . As 
our data shows, the true asymptopic value of V is reached 
only at much larger values oi L y . At T = 0.7 J /ks arid 
/ = 0.1/ c , for example, the asymptopic value is larger 
than V m by five-fold. 

In the remaining part of the paper, we propose an ex- 
planation of the unusual finite-size effect, based on the 
standard picture of creation and annhilation of free vor- 
tices under an applied current / For small /, the 
supercurrent through the network is maintained by an 
average phase jump A(f> = I / I c per column. Passage 
of a free vortex across the strip (i.e., in the y-direction) 
decreases the overall phase advance along the strip by 
an amount 2ir. In a steady-state situation, this loss is 
compensated by a difference in the phase velocities at 
the two opposing ends of the strip, and hence a voltage 
drop across the system. Simple counting then yields the 
following expression for the average voltage drop per col- 
umn in a strip of length L x 1 (see, e.g. Ref. pj), 

V = j^-(N + v + -JV_tJ_). (5) 

Here N± are the number of ± vortices in the system, and 
v± are their respective mean velocities in the y-direction. 
Although Eq. (§) does not discriminate between bound 
and free vortices, it is clear that correlated drift of a 
vortex-antivortex pair does not contribute to V. 

In the AHNS theory, the energy of a vortex/antivortex 
pair of size R and oriented perpendicular to the external 
current is given by, 

E(R) = 2E c + 2Ttj(\nR- jR), (6) 

where E c is the core energy of a single vortex. Equation 
@ has a maximum E = Eb = 2E C + 2ir J\n(I c /I) — 2-kJ 
at R = ?b = Ic/I- This is the activation energy for pair 
breaking. The rate for such processes follows the Arrhc- 
nius law, 

r ~ e - E B/k B T _ (J// C )2X. (7) 

In the second step of the argument, one writes down a 
mean-field rate equation governing the density p — p± of 
free vortices, 

p = T- 1P \ (8) 



where 7 is the diffusion constant. The steady-state vor- 
tex density is then given by p^ = (r/7) 1 / 2 . Going back 
to Eq. (0), one finds 

Voo ~PocV± ~I X+ \ (9) 

where we have assumed the drift velocity v± of a free 
vortex to be proportional to I. 

In the RSJ dynamics, the spin-waves have a 
wavelength-independent relaxation time, of the order of 
the basic time unit afi/2el c . Hence the equilibrium in- 
teraction between vortices in (^|) is un-retarded in the 
dynamic case. Equation (||), on the other hand, is based 
on the assumption that there is no spatial correlation 
between vortices and antivortices. It turns out that 
this assumption is violated when the width of the strip 

— 1/2 

is less than a certain characteristic size l r — p^ ~ 
(I c /I) x / 2 . In a bulk system, there is typically one free 
vortex/antivortex pair in an area of linear size l r at any 
given time. This is no longer the case when L y < l r . 
Under periodic boundary conditions in the y-direction, 
a vortex may annhilate with an antivortex created from 
the same bound pair. Consequently, the life time of the 
two are shortened, giving rise to a spurious lower voltage 
as seen in Fig. 2. 

To obtain a semi-quantitative estimate of the voltage 
reduction, let us consider a simplified model of vortex- 
antivortex recombination. In a coarse-grained descrip- 
tion, we take lb to be the basic unit of length. The 
time for a single vortex to traverse the width of the 
strip is t = L y /v± = L y l\,. During this period, the 
typical displacement of the vortex in the x-direction is 
I = t 1 / 2 . The probability for the vortex to be within 
a distance from its initial position after one round is 
thus pi = = (lb/ ' Ly) 1 / 2 . This is the probability for a 
vortex-antivortex pair to recombine in one round across 
the strip. Assuming pi C 1, the probability for this to 
happen in two rounds is P2 — (lb/2L y ) 1 / 2 . In general, 
p n ~ pin -1 / 2 provided pin 1 / 2 < 1. The survival proba- 
bility after n rounds is 1 — X)"=iP« — 1 — pin 1 / 2 . This 
yields a self-recombination time, 

T s ^p^ 2 T = L 2 y . (10) 

The density of free vortices in this case is p s ~ Tt s ~ 
L 2 y (I/I c ) 2 *. From Eq. (g) we obtain, 

V s ~ p s v ± ~ L 2 y (I/I c fx+\ (11) 

which increases with increasing L y . 

The typical displacement of the vortex in the x- 
1/2 

direction over t s is r s = L y . To preempt self- 
recombination, we need TL 2 t s > 1, i.e., creation of the 
second vortex/antivortex pair in the same region in space 
within the time r s . This yields L y = (I c /I) x ^ 2 = l T \ 
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Therefore Eq. ( Jll| ) applies in the intermediate regime 
l h < L y < / r . 

The situation is somewhat different when L a < lb- 
In this case, the pair-breaking interaction in (^) due to 
the external current is always weaker than the vortex- 
antivortex attraction and hence can be treated as a per- 
turbation. The I-V curves can in principle be calculated 
in a linear response theory [Q. Qualitatively, though, 
we may estimate the finite-size effect from the equilib- 
rium relaxation time of a vortex-antivortex pair on scale 



' eq 



(12) 



where z = 2\ — 2 is the dynamic exponent. We now di- 
vide the strip into square segments of L y columns each. 
The rate of creating a vortex-antivortex pair of size L y in 
such a segment is . In thermal equilibrium, the vor- 
tex/antivortex may drift in either direction around the 
strip and recombine to generate a ±2n phase jump. This 
symmetry is broken by the external current which en- 
hances the —2n jump (along the current direction) by a 
factor equal to the ratio between the energy difference of 
the two possibilities and fc^T. This leads to the estimate, 



Vohr 



1 



IrknT 



{I/Io)L-- 



(13) 



In the Ohmic regime, V decreases rapidly with increasing 
L y , as seen in Fig. 2. 

At L y — lh, the two expressions (JTTJ) and (|l^) coincide 
to give the minimum value V m = (I/I c ) 1+Z . This is in- 
deed the power-law form obtained by Minnhagen et al. 
by assuming lb as the only relevant length scale. 

A direct consequence of the self-recombination picture 
proposed above is the intermittent phase difference jumps 
across a square array of size less than l r . In the steady 
state, the phase difference A<fi across the system, aver- 
aged in the direction perpendicular to the current, is ex- 
pected to grow linearly with time. However, if there is 
only at most one pair of free vortex/antivortex present 
in the system at any given time, the growth of A(j) is in- 
termittent. This is indeed the case for L y < l r , where the 
life time of a free vortex/antivortex pair is less than the 
time it takes to create the pair. The step-wise growth 
of A(j) has been observed in a previous numerical study 
using a square array J^] which we have confirmed. A sim- 
ilar phenomenon occurs in the L y < lb regime, although 
in this case there are occasionally reverse steps. For a 
long strip, however, A<f> is a much smoother function of 
time due to averaging along the strip. 

In summary, extensive simulations of the Joscphson 
junction array in strip geometry reveal the existence of 
three size regimes with distinct current- voltage relation- 
ships. Arrays of size less than lb ~ I^ 1 are in the 
Ohmic regime and obey dynamic scaling as proposed 
by Minnhagen et al. Between lb and a second scale 



l T ~ /~ x / 2 , a reversed size-effect is observed. The power- 
law I-V dependence of AHNS is shown to hold in the 
asymptotic regime. Our study reconciles a long-standing 
dispute between two theoretical approaches, and provides 
a framework for systematic analysis of finite-size effects 
in 2D superfluid and superconducting systems with equal 
number of vortices and antivortices. 
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